Prediction of the bottomonium D-wave spectrum from full lattice QCD 
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We calculate the full spectrum of D-wave states in the T system in lattice QCD for the first time, 
using an improved version of NonRelativistic QCD on coarse and fine 'second generation' gluon field 
configurations from the MILC collaboration that include the effect of up, down, strange and charm 
quarks in the sea. Taking the 2S — IS splitting to set the lattice spacing, we determine the 3 D2 — IS 
splitting to 2.3%, and find agreement with experiment. Our prediction of the fine structure relative 
to the 3 D 2 gives the 3 D 3 at 10.181(5) GeV and the 3 Di at 10.147(6) GeV. We also discuss the 
overlap of 3 Di operators with 3 Si states. 



Introduction. The spectrum of bb states has provided 
a very important testing ground for strong interaction 
physics because of the number of radial and orbital exci- 
tations that are 'gold-plated', i.e. well below the thresh- 
old for decay to B mesons. The recent discovery of the 
r] b {lS) 1\ and h b (lP) and h b (2P) mesons [2] filled in 
important gaps in the spin-singlet states. The mass of 
the r\b meson had previously been predicted by lattice 
QCD [3] and the h b meson masses were widely expected, 
and found, to be very close to the spin-average of their 
associated spin-triplet states. 

The key missing gold-plated mesons are now the 
Y(l-D) states. These are very difficult to find experimen- 
tally although the 3 D 2 has been seen in radiative decay 
from the T(3S') [I]. Masses of the D-wave states have 
been predicted in potential model calculations (see, for 
example |6] ) , but it is hard to quantify the errors in 
these predictions except by using different forms for the 
potentials. 

In lattice QCD the starting point is QCD itself. The 
difficulties with the D-wave states then stem from the sig- 
nal to noise ratio; the signal falls exponentially in lattice 
time with the D-wave mass but the noise falls with the 
smaller ground state S-wave mass. Very large samples 
of meson correlators on full QCD gluon field configura- 
tions are then needed to obtain a reliable signal. Here 
we give the first results from lattice QCD that are able 
to distinguish the fine structure of D-wave states. 

Lattice Calculation. We use 'second generation' gluon 
field configurations recently generated by the MILC col- 
laboration [7j . These have a gluon action fully improved 
through a s a 2 [5] and include the effect of u, d, s and 
c quarks in the sea using the Highly Improved Stag- 
gered Quark formalism [D]. The u and d quarks have 
the same mass, mi, so the configurations are denoted as 
rif = 2 + 1 + 1. We use three ensembles to give two 
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values of the lattice spacing and two values of mi. The 
parameters of the ensembles are given in Table [TJ we la- 
bel them as 3, 4 and 5 from earlier work [lOj in which we 
mapped out the S and P-wave bottomonium spectrum 
and determined the lattice spacing from the T (25 — IS) 
splitting. 

We calculate b quark propagators on these configura- 
tions using an improved lattice discretisation of NonRel- 
ativistic QCD (NRQCD). NRQCD is an expansion in 
powers of the heavy quark velocity and therefore good 
for b quarks since v 2 /c 2 « 0.1 inside their bound states. 
The Hamiltonian includes all terms through 0(v 4 ) [ID] : 
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Here V is the symmetric lattice derivative and is the 
lattice discretization of the continuum D™. E and B 
are the chromoelectric and chromomagnetic fields, omj, 
is the bare b quark mass, which is tuned by determination 
on the lattice of the spin-average of ground-state T and 
r/ b meson masses. This was done in |10j to give the values 
used here, quoted in Table [Tl] 

The i> 4 terms in SH have coefficients Ci whose values 
are fixed from matching lattice NRQCD to full QCD, ei- 
ther perturbatively or nonperturbatively. Here we use 
coefficients for cj., C5 and c§ that include 0(a s ) correc- 
tions, as described in [lfj . The coefficients C3 and C4 
of the spin-dependent i> 4 terms have been tuned from a 
study of the fine structure of the Xfe(l-P) states. We find 
C3 = 1.0 with an error of 0.1. C4 is significantly larger. 
Here we use C4 = 1.25 on the coarse lattices and 1.10 on 
the fine lattices. These agree within 0.1 both with the 
value required to give P-wave fine structure in agreement 
with experiment and with the 0(a s ) improved result [10) . 

To make meson correlators for D-wave states we use a 
quark propagator made from either a local or a smeared 



2 



TABLE I. Details of the MILC gluon field ensembles used in 
this paper, a is the lattice spacing in fm determined from the 
T (25* — IS) splitting and L/a and T/a give the lattice size. 
ami, arris and am c are the sea quark masses in lattice units. 
Ensembles 3 and 4 are denoted "coarse" and 5, "fine." 



Set 


a (fm) 


ami 


am s 


am c 


L/a x T/a 


3 


0.1219(9) 


0.0102 


0.0509 


0.635 


24x64 


4 


0.1195(10) 


0.00507 


0.0507 


0.628 


32x64 


5 


0.0884(6) 


0.0074 


0.037 


0.440 


32x96 



source which has appropriate derivatives applied to it to 
generate a _D-'wavefunction'. This propagator is then 
combined with a local propagator and the same deriva- 
tives and smearings applied at the sink to create a 2 x 2 
matrix of correlators for each D-w&ve state. The com- 
plete set of combinations of spin matrices and derivatives 
needed is given in [II] . Note that the spin-2 and spin- 
3 representations split into irreducible representations 
of the lattice rotational group {Ai, A 2 , E, Ti, T 2 }, which 
must be considered independently since their masses can 
differ by discretisation errors. Very high statistics is re- 
quired - we have typically 32,000 correlators for every 
source operator per ensemble, using multiple time sources 
per configuration. The time sources are binned over for 
analysis. 

Bayesian fitting [T^] is used to extract the spectrum 
from the correlators using fit function: 



Gmcson(n sc ,n sk ;t) = ^2a(n sc ,k)a*(n sk ,k)e Ekt . (2) 

E k is the energy of the (k — l)th radial excitation and 
a(n, k) label the amplitudes depending on source and sink 
smearing. We fit all the D-wave states together taking 
the 3 D 2 E state as the reference state, with a prior of 
width 0.1 on its ground-state energy. Relative to that we 
take prior value ± 40 MeV on the ground-state energy 
of the other states. We take priors 0.5 ± 0.5 GeV on 
radial excitation energies and 0.1 ± 1.0 on amplitudes. 
We fit correlators from time t /a = 2 to 12 except for the 
local-local correlators which we take from t/a — 9 to 12. 

Results. The results from our fits for each D-wsve 
lattice representation on each ensemble are given in Ta- 
ble [II] We use 7i cxp = 3 on sets 3 and 4 and n cxp = 4 
on set 5 since these have the highest posterior prob- 
ability [12]; values and errors have stabilised at this 
point and x 2 /dof < 1. We also give the ratio Rd — 
(1 3 L> 2 - lSV(2 3 Si - l 3 Si) where IS is the spin-average 
of T and rjb energies from [10] and 1 3 D2 is the dimension- 
weighted average of the lattice 3 D 2 e and 3 D 2 t 2 results. 

Rd is plotted along with similarly defined Rs and Rp 
from [TU] in Figure [T] To obtain a physical result for Rd 
we fit to the same form used in |10j for Rs and Rp, allow- 
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FIG. 1. Results for the ratio of the 1 3 Z?2 — IS splitting to the 
2S — IS splitting in the T system plotted against the square 
of the lattice spacing determined from the 2S — IS splitting. 
Other ratios from [10] are shown for comparison. The grey 
shaded bands give the physical result obtained from a fit to 
the data as described in the text, and with the full error of 
Table [Till The black open circles slightly offset from a = are 
from experiment |13| . 

ing for lattice spacing and sea quark mass dependence: 
R = R phys [l + 2hSxi(l + Cl (aAf) (3) 

+ X] c i( aA ) 2 ' J ( 1 + CjbSx m + Cjbb(Sx m ) 2 )]- 

i=i,2 

Here Sxi is (am/ / am s ) — (rni/m s ) p hys for each ensemble. 
(m;/m s )ph ys is taken from lattice QCD as 27.2(3) [14] . 
5x m — (anib — 2.65)/1.5 allows for arrib effects from 
NRQCD in the discretisation errors over our range of 
a values. A, taken as 500 MeV, sets the scale for phys- 
ical a-dependence. Fit priors are as in [10]: 1.0(0.5) on 
-Rphys! 0.0(0.3) on a 2 terms; 0.0(1.0) on higher order in a; 
0.0(0.015) on bi. The physical result we obtain for Rd is 
1.307(30), after adding an additional NRQCD systematic 
error for missing v e terms |10j . This is to be compared 
to the experimental value of 1.280(3). A complete error 
budget for Rd is given in Table jlLT) 

In Figure [2] we plot the masses of all the lattice repre- 
sentations relative to the spin average of all 1 3 D states 
for coarse set 4, using the 2S — 15 splitting to set the 
scale (Table [l| . We see that the lattice representations 
for each spin agree well with each other within our size- 
able statistical errors. The hyperfine splitting, between 
the 1 D 2 and the spin average of 3 D states is expected 
to be very small, following results for P-wave states. We 
find it to be zero to within 10 MeV. 

Figure [3] shows the results from all three sets, using a 
dimension- weighted average of results, including the cor- 
relations from the fit, for the different lattice represen- 
tations for the 3 £>3 and ^ 1,3 'l)2- Results are consistent 
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TABLE II. Fitted D-wave energies on each ensemble. Errors 
are from statistics/fitting only. C3 = 1.0 on all ensembles, 
C4 = 1.25 on sets 3 and 4 and 1.10 on set 5. aA(x) — aE(x) — 
aE( 3 D). Rx and Ax are defined in the text. The A 2 irrep. 
on set 5 is fit separately and not included in the splittings. 



Set 3 

amb = 2.66 



Set 4 

amb = 2.62 



Set 5 

amb = 1.91 



aE{l L D 2E ) 


0.705(10) 


0.694(12) 


0.594(5) 


aS(l 1 D 2 T 2 ) 
aE(l 3 D 1Tl ) 


0.711(8) 


0.693(10) 


0.589(3) 


0.695(7) 


0.680(10) 


0.575(8) 


aE{l 3 D 2E ) 


0.698(10) 


0.692(10) 


0.588(4) 


aE(l 3 D 2T2 ) 


0.702(8) 


0.691(10) 


0.589(4) 


aE(l 3 D 3A2 ) 
aE(l 3 D 3Tl ) 


0.707(10) 


0.704(10) 


0.597(4) 


0.715(7) 


0.705(8) 


0.596(4) 


aE(l 3 D 3T2 ) 


0.714(7) 


0.696(9) 


0.594(3) 



aA(l 1 D 2 ) 


0.0029(31) 


0.0004(37) 


0.0027(27) 


aA(l 3 Di) 


-0.0104(34) 


-0.0137(44) 


-0.0137(62) 


aA(l 3 D 2 ) 


-0.0047(23) 


-0.0021(21) 


0.0001(20) 


aA(l 3 D 3 ) 


0.0078(22) 


0.0074(27) 


0.0069(20) 


Rd 


1.318(23) 


1.303(26) 


1.309(16) 


aA L .s 


0.0038(11) 


0.0040(13) 


0.0037(13) 




-0.0005(9) 


0.0009(9) 


0.0016(15) 


Rl-s 


0.44(13) 


0.49(17) 


0.60(21) 


R Sij 


-0.26(52) 


0.53(50) 


1.1(1.0) 


TABLE III. 


Complete error 


budget for Rd 


in %. Finite 


volume and m b tuning errors 


are negligible. 




Rd 



stats/fitting 

a-dependence 

mj-dependence 

NRQCD am^-dependence 

NRQCD systematics 

electromagnetism / r\b annihilation 

Total 



1.4 
1.4 
0.5 
0.1 
1.0 
0.2 
2.3% 




FIG. 2. Results for the separate irreducible representations of 
the lattice rotation group making up each continuum D-wave 
state on coarse set 4. 



30 



20 



I 



Result 1 
Set 3 1 — fc- 
Set 4 1— -B- 
Set5 ' h- 



FIG. 3. D-wave masses plotted relative to the D spin-average 
for all sets using the 2S — IS splitting to set the scale. The 
red shaded bands show our final results using ratios of com- 
binations of splittings to those of P-wave states as described 
in the text. 



between the fine and coarse sets and between different 
sea light quark masses for the two coarse sets. 

To arrive at a final result for D-wave fine structure we 
study combinations of 3 D spin-splittings that are sensi- 
tive either to an L • S or to a tensor Sy interaction (S • S 
takes the same value for all 3 D states). Writing 

Mj = M( 3 D) + Al s {L ■ S) + A* (S tj ) (4) 

gives 

Af. s = (14M 3 - 5M 2 - 9M x )/60 

Ag. = -7(2M 3 - 5M 2 + 3M0/120. (5) 

Table [TT] gives our results for these splittings. In Fig- 
ure [J] we plot ratios to the equivalent 3 P splitting com- 
binations: R x = A x /A x with A£. s = (5M 2 - 3Mi - 
2M )/12 and Af = -5(M 2 - 3M X + 2M )/72. Val- 
ues for A p for these ensembles are given in [10] (without 
factors of 1/12 and -5/72). The experimental values are 
A£<j = 13.65(27) MeV and Af.. = 3.29(9) MeV [T3]. 



The advantage of using these combinations is that they 
depend purely on one of the spin-dependent coefficients 
of the NRQCD action. On set 5 we did not use exactly 
the same values for C3 and C4 in our study of P and D 
waves. However we can correct for this in Figure [4] since 
As w ex C4 and A^.g oc C3. Once this slight adjustment 
is done the dependence on C3 4 cancels between P and D 
states and so errors from the uncertainty in these coeffi- 
cients are much reduced. 

We fit the fine-structure R values to the same form 
used earlier in eq. [4] to extract physical results: 

R L . S = 0.49(11); R SlJ . = 0.26(35). (6) 

We have included an additional systematic error of 10% 
to allow for missing v e terms from our NRQCD action but 
the lattice statistical error dominates. We then combine 
the R values with experimental results from IP levels to 
give the following 3 D splittings: 

3 D 3 - 3 Di = 34(8)MeV 
3 D 3 - 3 D 2 = 18(5)MeV 
3 D 2 - 3 Di = 17(6)MeV. (7) 
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FIG. 4. Ratios iisy and 7?l s (multiplied by -5 for clarity) 
plotted against the square of the lattice spacing. The grey 
bands give our physical results. 
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FIG. 5. Correlators made from different combinations of local 
3 Si and :i Di operators at source and sink, plotted with a 
logarithmic y axis as a function of lattice time, t. Results are 
from set 3. 



Our fine structure splittings are somewhat larger than 



typical results from potential models [5] [BJ, where the 
3 D 3 - 3 L»i splitting lies in the range 10-20 MeV. This 
can be traced to a larger value for Rl-s than is obtained, 
for example, in [5], based on specific forms for the spin- 
dependent potentials. 

One issue that we have neglected above is that the 3 Di 
state has J PC = 1~~ in common with 3 5i states. On 
the lattice, in principle, any operator with 1 quantum 
numbers will be able to create all 1 states. In practice 
the amplitude for 3 Si states to be created by the opera- 
tors that we use for the 3 Di is very small and vice versa. 
We illustrate that in Fig. [5] where we show correlators 
from set 3 that use a local iSi or 3 £>i operator at source 
and sink compared to the cross-correlator that has a lo- 
cal 3 Si operator at the source and 3 -Di at sink or vice 
versa. The cross-correlator is much smaller in magnitude 
than either of the diagonal correlators at small t values. 
The exponential fall-off (as seen in the slope of the log 
plot) of the cross-correlator matches that of the 3 S\ cor- 
relator at large times, where the 3 Z?i correlator fall-off is 
dominated by that of the heavier 3 Di state. If we fit the 
complete set of 3 Si and 3 D\ correlators together, includ- 
ing the local cross-correlators of Fig. [5j we obtain results 
in agreement with our separate fits for 3 S\ (in [TUj) and 
3 Di masses. We also find, for example, that the ampli- 
tude a( 3 Dij oca i, T) from eq. [2]is 0.0052(1) times that of 
a( Sij oca i, T). 

Conclusions. We give the first full lattice QCD results 
for the D-wave states of bottomonium including the fine 
structure. We obtain a mass of 10.179(17) GeV for the 
1 3 _D 2 to be compared with 10.1637(14) GeV from exper- 
iment [13] . Using the experimental result for the 1 3 L>2 
mass we predict masses of 10.181(5) GeV for the 1 3 D 3 , 
10.147(6) GeV for the l 3 L>i and 10.169(10) for the 1 D 2 . 
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